/*
 * File:   main.cpp
 * Author: thiago
 *
 * Created on 8 de Novembro de 2009, 22:26
 */

//#define DEBUG true
//#define DEBUG_VALORES true
//#define DEBUG_OPENGL false

#include <stdlib.h>
#include <malloc.h>
#include <cmath>
#include <iostream>
#include <math.h>

//Includes GL
#include <GL/glew.h>
#include <GL/glut.h>
#include <GL/glext.h>

//Includes CUDA
#include <cuda_runtime.h>
#include <cutil.h>
#include <vector_types.h>
#include <cuda_gl_interop.h>
#include <cutil_gl_error.h>
#include <assert.h>
#include <cudpp/cudpp.h>

#include "Particula.h"

#include "RadixSort/radixsort.h"
#include "Scan/Scan.h"

#define ABS(a)		(((a) < 0) ? -(a) : (a))

using namespace std;

extern "C" void launch_CalculaIDGridCell(unsigned int blocksPerGrid, unsigned int threadsPerBlock,
        float3* posicoes, unsigned int* values, unsigned int* keys, int numeroParticulas,
        int numColunas, float deltaX, float deltaY, unsigned int* count, float XminGrid, float YminGrid);

extern "C" void launch_FindNeighbors(unsigned int blocksPerGrid, unsigned int threadsPerBlock,
        float3* posicoes, Particula* particulas, Particula* particulas2, unsigned int* values, unsigned int* keys,
        unsigned int* count, unsigned int* offSet, int numeroColunasGridCell, int numeroLinhasGridCell, int numeroParticulas,
        float deltaX, float deltaY, float XminGrid, float YminGrid, float massaEspecificaInicial, float massa);

extern "C" void launch_ZeraCountEOffSet(unsigned int blocksPerGrid, unsigned int threadsPerBlock, int numeroCelulas,
        unsigned int* count, unsigned int* offSet);

float massa;
float massaEspecificaInicial;
const float constGas = 0.8;
const float coeficienteViscosidade = 0.00005;
float passoIntegracao = 0.8;
//Vetores para controle das particulas.
float3* posicoesParticulas;
float3* posicoesParticulas2;
float3* posicoesParticulasT;
float3* posicoesParticulasT2;
float3* forcaPressaoVetor;
float3* forcaAceleracaoVetor;
float3* forcaViscosidadeVetor;
Particula* particulas;
Particula* particulasDevice;
Particula* particulas2;
Particula* particulasDevice2;
Particula* particulasT;
Particula* particulasT2;
unsigned int* keys;
unsigned int* keysDevice;
unsigned int* values;
unsigned int* valuesDevice;
unsigned int* count;
unsigned int* offSet;
unsigned int* countDevice;
unsigned int* offSetDevice;
int numeroParticulas;
int numeroColunasGridCell;
int numeroLinhasGridCell;
int numeroCelulasGrid;
float XminGrid = 0.0, XmaxGrid = 1.0, YminGrid = 0.0, YmaxGrid = 1.0, deltaXGrid, deltaYGrid;
float recuoXGrid = 0.02, recuoYGrid = 0.02;
//float* deltaXGridDevice;
//float* deltaYGridDevice;
int threadsPerBlock, blocksPerGrid;
//Utilizado no desenho das particulas.
int colunas, linhas;
float posicaoX, posicaoY, passoX, passoY;
bool mostraGrid, mostraEixo, mostraVetores, mostraPressao, mostraViscosidade, mostraAceleracao;
//VBO
GLuint positionsVBO;
float3* VBO;
GLuint positionsVBO2;
float3* VBO2;
bool VBO_Ativo;
const bool gpu = false;

/*******************************************************************************
 *
 * METODOS DE DESENHO
 *
 ******************************************************************************/

void desenhaEixos() {

    glBegin(GL_LINES);
    glColor3f(1.0, 0.0, 0.0);
    glVertex3f(1.0, 0.002, 0.0);
    glVertex3f(0.0, 0.002, 0.0);
    glVertex3f(1.0, 0.001, 0.0);
    glVertex3f(0.0, 0.001, 0.0);

    glColor3f(0.0, 1.0, 0.0);
    glVertex3f(0.001, 1.0, 0.0);
    glVertex3f(0.001, 0.0, 0.0);
    glVertex3f(0.000, 1.0, 0.0);
    glVertex3f(0.000, 0.0, 0.0);

    //		glColor3f(0.0, 0.0, 1.0);
    //		glVertex3f(0.0, 0.0, 1.0);
    //		glVertex3f(0.0, 0.0, 0.0);

    glEnd();

}

void desenhaGrid() {

    float posicaoI = YmaxGrid, posicaoJ = XminGrid;

    glBegin(GL_LINES);

    glColor3f(1.0, 0.4, 0.0);
    for (int i = 0; i <= numeroLinhasGridCell; i++) {
        glVertex3f(XminGrid, posicaoI, 0.0);
        glVertex3f(XmaxGrid, posicaoI, 0.0);
        posicaoI = posicaoI - deltaYGrid;
    }
    for (int i = 0; i <= numeroColunasGridCell; i++) {
        glVertex3f(posicaoJ, YmaxGrid, 0.0);
        glVertex3f(posicaoJ, YminGrid, 0.0);
        posicaoJ = posicaoJ + deltaXGrid;
    }

    glEnd();
}

void desenhaVetoresAceleracao(){
    glBegin(GL_LINES);

    glColor3f(1.0, 0.0, 0.0);
    for (int i = 0; i <= numeroParticulas; i++) {
        float Ex, Ey, PontoX, PontoY;
        if (forcaAceleracaoVetor[i].x != 0) {
            Ex = passoX / (5 * forcaAceleracaoVetor[i].x);
        } else {
            Ex = 0;
        }
        if (forcaAceleracaoVetor[i].y != 0) {
            Ey = passoY / (5 * forcaAceleracaoVetor[i].y);
        } else {
            Ey = 0;
        }
        PontoX = ABS(Ex) * forcaAceleracaoVetor[i].x;
        PontoY = ABS(Ey) * forcaAceleracaoVetor[i].y;
        glVertex3f(posicoesParticulasT[i].x, posicoesParticulasT[i].y, 0.0);
        glVertex3f(posicoesParticulasT[i].x + PontoX,
                posicoesParticulasT[i].y + PontoY,
                0.0);
    }
    glEnd();

}

void desenhaVetores(){
    glBegin(GL_LINES);

    glColor3f(1.0, 1.0, 0.0);
    for (int i = 0; i <= numeroParticulas; i++) {
        float diferencaX = (posicoesParticulasT2[i].x - posicoesParticulasT[i].x);
        float diferencaY = (posicoesParticulasT2[i].y - posicoesParticulasT[i].y);
        float Ex, Ey, PontoX, PontoY;
        if (diferencaX != 0) {
            Ex = passoX / (5 * diferencaX);
        } else {
            Ex = 0;
        }
        if (diferencaY != 0) {
            Ey = passoY / (5 * diferencaY);
        } else {
            Ey = 0;
        }
        PontoX = ABS(Ex) * diferencaX;
        PontoY = ABS(Ey) * diferencaY;
        glVertex3f(posicoesParticulasT[i].x, posicoesParticulasT[i].y, 0.0);
        glVertex3f(posicoesParticulasT[i].x + PontoX,
                posicoesParticulasT[i].y + PontoY,
                0.0);
    }

    glEnd();
}

void desenhaVetoresPressao(){
    glBegin(GL_LINES);

    glColor3f(1.0, 1.0, 1.0);
    for (int i = 0; i <= numeroParticulas; i++) {
        float Ex, Ey, PontoX, PontoY;
        if (forcaPressaoVetor[i].x != 0) {
            Ex = passoX / (5 * forcaPressaoVetor[i].x);
        } else {
            Ex = 0;
        }
        if (forcaPressaoVetor[i].y != 0) {
            Ey = passoY / (5 * forcaPressaoVetor[i].y);
        } else {
            Ey = 0;
        }
        PontoX = ABS(Ex) * forcaPressaoVetor[i].x;
        PontoY = ABS(Ey) * forcaPressaoVetor[i].y;
        glVertex3f(posicoesParticulasT[i].x, posicoesParticulasT[i].y, 0.0);
        glVertex3f(posicoesParticulasT[i].x + PontoX,
                posicoesParticulasT[i].y + PontoY,
                0.0);
    }

    glEnd();
}

void desenhaVetoresViscosidade(){
    glBegin(GL_LINES);

    glColor3f(1.0, 0.6, 0.0);
    for (int i = 0; i <= numeroParticulas; i++) {
        glVertex3f(posicoesParticulasT[i].x, posicoesParticulasT[i].y, 0.0);
        glVertex3f(posicoesParticulasT[i].x + 100000*forcaViscosidadeVetor[i].x,
                   posicoesParticulasT[i].y + 100000*forcaViscosidadeVetor[i].y,
                   0.0);
    }

    glEnd();
}

//Constroi o VBO
void buildVBO() {

#ifdef DEBUG_GL
    cout << "CONSTRUINDO VBO 1!" << endl;
    cout << "ANTES DO GenBuffer!" << endl;
    glGenBuffers(1, &positionsVBO);
    cout << "ANTES DO BindBuffer com o positionVBO!" << endl;
    glBindBuffer(GL_ARRAY_BUFFER, positionsVBO);
    cout << "ANTES DO BufferData!" << endl;
    unsigned int size = numeroParticulas * sizeof (float3);
    glBufferData(GL_ARRAY_BUFFER, size, posicoesParticulas, GL_DYNAMIC_DRAW);
    cout << "ANTES DO BindBuffer com o 0!" << endl;
    glBindBuffer(GL_ARRAY_BUFFER, 0);

    cout << "CONSTRUINDO VBO 2!" << endl;
    cout << "ANTES DO GenBuffer2!" << endl;
    glGenBuffers(1, &positionsVBO2);
    cout << "ANTES DO BindBuffer com o positionVBO2!" << endl;
    glBindBuffer(GL_ARRAY_BUFFER, positionsVBO2);
    cout << "ANTES DO BufferData2!" << endl;
    unsigned int size2 = numeroParticulas * sizeof (float3);
    glBufferData(GL_ARRAY_BUFFER, size2, posicoesParticulas2, GL_DYNAMIC_DRAW);
    cout << "ANTES DO BindBuffer com o 0!" << endl;
    glBindBuffer(GL_ARRAY_BUFFER, 0);
#else
    glGenBuffers(1, &positionsVBO);
    glBindBuffer(GL_ARRAY_BUFFER, positionsVBO);
    unsigned int size = numeroParticulas * sizeof (float3);
    glBufferData(GL_ARRAY_BUFFER, size, posicoesParticulas, GL_DYNAMIC_DRAW);
    glBindBuffer(GL_ARRAY_BUFFER, 0);

    glGenBuffers(1, &positionsVBO2);
    glBindBuffer(GL_ARRAY_BUFFER, positionsVBO2);
    unsigned int size2 = numeroParticulas * sizeof (float3);
    glBufferData(GL_ARRAY_BUFFER, size2, posicoesParticulas2, GL_DYNAMIC_DRAW);
    glBindBuffer(GL_ARRAY_BUFFER, 0);
#endif

    CUDA_SAFE_CALL(cudaGLRegisterBufferObject(positionsVBO));

    CUDA_SAFE_CALL(cudaGLRegisterBufferObject(positionsVBO2));

    CUT_CHECK_ERROR_GL(); //cout << "glGetError(): " << glGetError() << endl;
}

void deleteVBO(){
    CUDA_SAFE_CALL(cudaGLUnregisterBufferObject(positionsVBO));
    CUDA_SAFE_CALL(cudaGLUnregisterBufferObject(positionsVBO2));
    glDeleteBuffers(1, &positionsVBO);
    glDeleteBuffers(1, &positionsVBO2);
}

// Função callback chamada para fazer o desenho
void display(void) {

    // Limpa a janela de visualização com a cor de fundo especificada
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    if (mostraEixo) {
        desenhaEixos();
    }

    if (mostraGrid) {
        desenhaGrid();
    }

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

#ifdef DEBUG_GL
    if (VBO_Ativo == true) {
        cout << "HABILITOU O GL_VERTEX_ARRAY" << endl;
        glEnableClientState(GL_VERTEX_ARRAY);
        // Especifica que a cor corrente é vermelha
        //         R     G     B
        glColor3f(0.0f, 0.6f, 1.0f);
        cout << "ANTES DO BindBuffer com o positionVBO!" << endl;
        glBindBuffer(GL_ARRAY_BUFFER, positionsVBO);
        cout << "ANTES DO GLVERTEXPOINTER" << endl;
        glVertexPointer(3, GL_FLOAT, 0, 0);
        cout << "ANTES DO GLDRAWARRAYS" << endl;
        glDrawArrays(GL_POINTS, 0, numeroParticulas);
        cout << "ANTES DO BindBuffer com o 0!" << endl;
        glBindBuffer(GL_ARRAY_BUFFER, 0);
        glDisableClientState(GL_VERTEX_ARRAY);
    } else {
        cout << "HABILITOU O GL_VERTEX_ARRAY" << endl;
        glEnableClientState(GL_VERTEX_ARRAY);
        // Especifica que a cor corrente é vermelha
        //         R     G     B
        glColor3f(0.0f, 0.6f, 1.0f);
        cout << "ANTES DO BindBuffer com o positionVBO!" << endl;
        glBindBuffer(GL_ARRAY_BUFFER, positionsVBO2);
        cout << "ANTES DO GLVERTEXPOINTER" << endl;
        glVertexPointer(3, GL_FLOAT, 0, 0);
        cout << "ANTES DO GLDRAWARRAYS" << endl;
        glDrawArrays(GL_POINTS, 0, numeroParticulas);
        cout << "ANTES DO BindBuffer com o 0!" << endl;
        glBindBuffer(GL_ARRAY_BUFFER, 0);
        glDisableClientState(GL_VERTEX_ARRAY);
    }
#else
    if (VBO_Ativo == true) {
        glEnableClientState(GL_VERTEX_ARRAY);
        // Especifica que a cor corrente é vermelha
        //         R     G     B
        glColor3f(0.0f, 0.6f, 1.0f);
        glBindBuffer(GL_ARRAY_BUFFER, positionsVBO);
        glVertexPointer(3, GL_FLOAT, 0, 0);
        glDrawArrays(GL_POINTS, 0, numeroParticulas);
        glBindBuffer(GL_ARRAY_BUFFER, 0);
        glDisableClientState(GL_VERTEX_ARRAY);
    } else {
        glEnableClientState(GL_VERTEX_ARRAY);
        // Especifica que a cor corrente é vermelha
        //         R     G     B
        glColor3f(0.0f, 0.6f, 1.0f);
        glBindBuffer(GL_ARRAY_BUFFER, positionsVBO2);
        glVertexPointer(3, GL_FLOAT, 0, 0);
        glDrawArrays(GL_POINTS, 0, numeroParticulas);
        glBindBuffer(GL_ARRAY_BUFFER, 0);
        glDisableClientState(GL_VERTEX_ARRAY);
    }
#endif

    if (mostraVetores) {
        desenhaVetores();
    }

    if (mostraPressao) {
        desenhaVetoresPressao();
    }
    if (mostraViscosidade) {
        desenhaVetoresViscosidade();
    }
    if(mostraAceleracao){
        desenhaVetoresAceleracao();
    }

    // Executa os comandos OpenGL
    glutSwapBuffers();

}

// Inicializa parâmetros de rendering
void init(void) {
    // Define a cor de fundo da janela de visualização como preta
    glClearColor(0.0f, 0.0f, 0.0f, 1.0f);
    glPointSize(3.0);
    /*
        glClearDepth(1.0);
        glEnable(GL_LIGHT0);
        glEnable(GL_LIGHTING);
        glEnable(GL_DEPTH_TEST);
     */

}

// Função callback chamada quando o tamanho da janela é alterado
void AlteraTamanhoJanela(GLsizei w, GLsizei h) {
    // Evita a divisao por zero
    if (h == 0) h = 1;

    // Especifica as dimensões da Viewport
    glViewport(0, 0, w, h);

    // Inicializa o sistema de coordenadas
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();

    // Estabelece a janela de seleção (left, right, bottom, top)
    if (w <= h)
        gluOrtho2D(0.0f, 1.0f, 0.0f, 1.0f/**h/w*/);
    else
        gluOrtho2D(0.0f, 1.0f/**w/h*/, 0.0f, 1.0f);
}

void printDeviceProps() {
    cudaDeviceProp prop;
    int dev;

    CUDA_SAFE_CALL(cudaGetDevice(&dev));
    cout << "dev = " << dev << endl;
    CUDA_SAFE_CALL(cudaGetDeviceProperties(&prop, dev));

    printf("Name:                    %s\n", prop.name);
    printf("Global Memory:           %d\n", prop.totalGlobalMem);
    printf("Shared Memory per Block: %d\n", prop.sharedMemPerBlock);
    printf("Register per Block:      %d\n", prop.regsPerBlock);
    printf("Warp Size:               %d\n", prop.warpSize);
    printf("Memory Pitch:            %d\n", prop.memPitch);
    printf("Max Threads Per Block:   %d\n", prop.maxThreadsPerBlock);
    printf("Max Threads Dimension:   %d %d %d\n", prop.maxThreadsDim[0], prop.maxThreadsDim[1], prop.maxThreadsDim[2]);
    printf("Max Grid Size:           %d %d %d\n", prop.maxGridSize[0], prop.maxGridSize[1], prop.maxGridSize[2]);
    printf("Constant Memory:         %d\n", prop.totalConstMem);
    printf("Version:                 %d.%d\n", prop.major, prop.minor);
    printf("Clock Rate:              %d\n", prop.clockRate);
    printf("Texture Alignment:       %d\n", prop.textureAlignment);
}

void printPosicoesParticulas() {
#ifdef DEBUG
    for (int j = 0; j < numeroParticulas; j++) {
        cout << "A posição (X) " << j << " do vetor A e igual a " << posicoesParticulas[j].x << endl;
    }
//    getchar();

    for (int j = 0; j < numeroParticulas; j++) {
        cout << "A posição (Y) " << j << " do vetor A e igual a " << posicoesParticulas[j].y << endl;
    }
//    getchar();
#endif
}

void printPosicoesParticulas2() {
#ifdef DEBUG
    for (int j = 0; j < numeroParticulas; j++) {
        cout << "A posição (X) " << j << " do vetor A e igual a " << posicoesParticulas2[j].x << endl;
    }
//    getchar();

    for (int j = 0; j < numeroParticulas; j++) {
        cout << "A posição (Y) " << j << " do vetor A e igual a " << posicoesParticulas2[j].y << endl;
    }
//    getchar();
#endif
}

void printDadosParticulas() {
#ifdef DEBUG
    for (int i = 0; i < numeroParticulas; i++) {
        cout << "id: " << values[i] << "\t idCell: " << keys[i] <<
                "\t Massa Esp: " << particulas[values[i]].massaEspecifica <<
                "\t Pressao: " << particulas[values[i]].pressao <<
                "\t Velo(x): " << particulas[values[i]].velocidade.x << endl;
    }
//    getchar();
#endif
}

void printDadosParticulas2() {
#ifdef DEBUG
    for (int i = 0; i < numeroParticulas; i++) {
        cout << "id: " << values[i] << "\t idCell: " << keys[i] <<
                "\t Massa Esp: " << particulas2[values[i]].massaEspecifica <<
                "\t Pressao: " << particulas2[values[i]].pressao <<
                "\t Velo(x): " << particulas2[values[i]].velocidade.x << endl;
    }
//    getchar();
#endif
}

void printCells() {
    for (int i = 0; i < numeroParticulas; i++) {
        cout << "id: " << values[i] << "\t idCell: " << keys[i] << endl;
    }
//    getchar();
}

void printCount() {
    for (int i = 0; i < numeroCelulasGrid; i++) {
        cout << "count[" << i << "] : " << count[i] << endl;
    }
//    getchar();
}

void radixSort() {
    int keybits = 32;

    // Creat the RadixSort object
    RadixSort radixsort(numeroParticulas, false);
    radixsort.sort(keysDevice, valuesDevice, numeroParticulas, keybits);

}

void scan() {
    Scan scanner(numeroParticulas);
    scanner.CriaOffset(countDevice, offSetDevice, numeroCelulasGrid);
}

void CalculaKernelsOnGPU() {
    if (VBO_Ativo) {
        // Map buffer object for writing from CUDA
        CUDA_SAFE_CALL(cudaGLMapBufferObject((void**) & VBO, positionsVBO));

        launch_ZeraCountEOffSet(blocksPerGrid, threadsPerBlock, numeroCelulasGrid, countDevice, offSetDevice);

        launch_CalculaIDGridCell(blocksPerGrid, threadsPerBlock, VBO, valuesDevice,
                keysDevice, numeroParticulas, numeroColunasGridCell, deltaXGrid, deltaYGrid, countDevice,
                XminGrid, YminGrid);

        scan();

        // Copy result from device memory to host memory
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        //cudaMemcpy(values, valuesDevice, numeroParticulas*sizeof(unsigned int), cudaMemcpyDeviceToHost);
        CUDA_SAFE_CALL(cudaMemcpy(count, countDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(offSet, offSetDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        printCells();

        printCount();

        cout << "\n Fazendo o radixSort" << endl;
#endif

        radixSort();
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(values, valuesDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        for (int i = 0; i < numeroParticulas; i++) {
            if (i < numeroCelulasGrid) {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << "\t offSet: " << offSet[i]
                        << "\t count[" << i << "] : " << count[i] << endl;
            } else {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << endl;
            }
        }
//        getchar();
#endif

        launch_FindNeighbors(blocksPerGrid, threadsPerBlock, VBO, particulasDevice, particulasDevice2, valuesDevice, keysDevice,
                countDevice, offSetDevice, numeroColunasGridCell, numeroLinhasGridCell, numeroParticulas,
                deltaXGrid, deltaYGrid, XminGrid, YminGrid, massaEspecificaInicial, massa);

        //CUDA_SAFE_CALL(cudaMemcpy(particulas, particulasDevice, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(particulas2, particulasDevice2, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        cout << "Imprimindo particulas t (particulas)." << endl;
        printDadosParticulas();
        cout << "\n Imprimindo particulas t + 1 (particulas2)." << endl;
        printDadosParticulas2();
#endif
        // Unmap buffer object
        CUDA_SAFE_CALL(cudaGLUnmapBufferObject(positionsVBO));

        VBO_Ativo = false;

    } else {

        // Map buffer object for writing from CUDA
        CUDA_SAFE_CALL(cudaGLMapBufferObject((void**) & VBO, positionsVBO));

        launch_ZeraCountEOffSet(blocksPerGrid, threadsPerBlock, numeroCelulasGrid, countDevice, offSetDevice);

        launch_CalculaIDGridCell(blocksPerGrid, threadsPerBlock, VBO, valuesDevice,
                keysDevice, numeroParticulas, numeroColunasGridCell, deltaXGrid, deltaYGrid, countDevice,
                XminGrid, YminGrid);

        scan();

        // Copy result from device memory to host memory
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        //cudaMemcpy(values, valuesDevice, numeroParticulas*sizeof(unsigned int), cudaMemcpyDeviceToHost);
        CUDA_SAFE_CALL(cudaMemcpy(count, countDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(offSet, offSetDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        printCells();

        printCount();

        cout << "\n Fazendo o radixSort" << endl;
#endif

        radixSort();
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(values, valuesDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        for (int i = 0; i < numeroParticulas; i++) {
            if (i < numeroCelulasGrid) {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << "\t offSet: " << offSet[i]
                        << "\t count[" << i << "] : " << count[i] << endl;
            } else {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << endl;
            }
        }
//        getchar();
#endif

        launch_FindNeighbors(blocksPerGrid, threadsPerBlock, VBO, particulasDevice2, particulasDevice, valuesDevice, keysDevice,
                countDevice, offSetDevice, numeroColunasGridCell, numeroLinhasGridCell, numeroParticulas,
                deltaXGrid, deltaYGrid, XminGrid, YminGrid, massaEspecificaInicial, massa);

        //CUDA_SAFE_CALL(cudaMemcpy(particulas, particulasDevice, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(particulas, particulasDevice, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));

#ifdef DEBUG
        cout << "Imprimindo particulas t + 1 (particulas)." << endl;
        printDadosParticulas();
        cout << "\n Imprimindo particulas t (particulas2)." << endl;
        printDadosParticulas2();
#endif
        // Unmap buffer object
        CUDA_SAFE_CALL(cudaGLUnmapBufferObject(positionsVBO));

        VBO_Ativo = true;
    }
}

float calculaMassaEspecifica(int i, int j, float h, int indiceAtual) {
    float massaEsp = 0;
    int idCell = (i * numeroColunasGridCell) + j;

        int kInicial = offSet[idCell];
        int kFinal = offSet[idCell] + count[idCell];

        //cout << "\t\t i  = " << i << "\t j = " << j << "\tidCell = " << idCell
        //       << "\t kInicial = " << kInicial << "\t kFinal = " << kFinal << endl;

        for (int k = kInicial; k < kFinal; k++) {
            int indiceParticulas = values[k];
            if (indiceParticulas != indiceAtual) {
                //cout << "\t\t Vizinho (massa esp.) => " << indiceParticulas << endl;
                float distancia = sqrt((posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x) * (posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x) + (posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y) * (posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y));
                //cout << "\t\tDistancia para o ponto " << indiceParticulas << " foi: " << distancia << "\t h = " << h << endl;
                if (distancia < h) {
                    //cout << "\t\t (Distancia menor que h)" << endl;
                    massaEsp += massa * (315 / (64 * 3.14 * pow(h, 9))) * pow((h * h - distancia * distancia), 3);
                    //cout << "\t\tmassa especifica da particula " << indiceParticulas << " eh " << massa * (315 / (64 * 3.14 * pow(h, 9))) * pow((h * h - distancia * distancia), 3) << endl;
                }
            }
        }

    //cout << "\t\t Valor da massa esp. foi " << massaEsp << endl;
    
    return massaEsp;
}

float3 calculaFocaPressao(int i, int j, float h, int indiceAtual) {
    float3 forcaPressao;
    forcaPressao.x = 0;
    forcaPressao.y = 0;
    forcaPressao.z = 0.0f;
    int idCell = (i * numeroColunasGridCell) + j;
    float h2 = h*h;

    int kInicial = offSet[idCell];
    int kFinal = offSet[idCell] + count[idCell];

    //cout << "\t\t i  = " << i << "\t j = " << j << "\tidCell = " << idCell
    //     << "\t kInicial = " << kInicial << "\t kFinal = " << kFinal << endl;


    for (int k = kInicial; k < kFinal; k++) {
        int indiceParticulas = values[k];
        if (indiceParticulas != indiceAtual) {
            //cout << "\t\t Vizinho (pressao) => " << indiceParticulas << endl;
            float distancia2 = (posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x) * (posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x) + (posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y) * (posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y);
            //float distanciaX = (posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x);
            //float distanciaY = (posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y);
            //float distanciaZ = (posicoesParticulasT[indiceAtual].z - posicoesParticulasT[indiceParticulas].z);
            assert(distancia2 >= 0);
            if (distancia2 < h2) {
                float distancia = sqrt(distancia2);
                //cout << "\t\t Distancia(" << distancia << ") menor que h(" << h << ")" << endl;
                float3 gradientePressao;
                float mod;
                mod = sqrt(pow(posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x, 2) +
                        pow(posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y, 2) +
                        pow(posicoesParticulasT[indiceAtual].z - posicoesParticulasT[indiceParticulas].z, 2));
                gradientePressao.x = (posicoesParticulasT[indiceParticulas].x - posicoesParticulasT[indiceAtual].x)/mod;
                gradientePressao.y = (posicoesParticulasT[indiceParticulas].y - posicoesParticulasT[indiceAtual].y)/mod;
                gradientePressao.z = (posicoesParticulasT[indiceParticulas].z - posicoesParticulasT[indiceAtual].z)/mod;

                forcaPressao.x += massa * ((particulasT[indiceParticulas].pressao + particulasT[indiceAtual].pressao) / 2 * particulasT[indiceParticulas].massaEspecifica) * gradientePressao.x * (15 / 3.14 * pow(h, 6) * (pow((h - distancia), 3)));
                //cout << "\t\t Pressao da particula " << indiceParticulas << " eh " << particulasT[indiceParticulas].pressao << "\t Massa esp. da particula " << indiceParticulas << " eh " << particulasT[indiceParticulas].massaEspecifica << endl;
                //cout << "\t\tforcaPressao.x da particula " << indiceParticulas << " eh " << forcaPressao.x << endl;
                forcaPressao.y += massa * ((particulasT[indiceParticulas].pressao + particulasT[indiceAtual].pressao) / 2 * particulasT[indiceParticulas].massaEspecifica) * gradientePressao.y * (15 / 3.14 * pow(h, 6) * (pow((h - distancia), 3)));
                //cout << "\t\tforcaPressao.y da particula " << indiceParticulas << " eh " << forcaPressao.y << endl;
                //forcaPressao.z += massa * ((particulasT[indiceParticulas].pressao + particulasT[indiceAtual].pressao) / 2 * particulasT[indiceParticulas].massaEspecifica) * (15 / 3.14 * pow(h, 6) * (pow((h - distanciaZ), 3)));
            }
        }
    }

    return forcaPressao;
}

float3 calculaFocaVisc(int i, int j, float h, int indiceAtual) {
    float3 forcaVisc;
    forcaVisc.x = 0;
    forcaVisc.y = 0;
    forcaVisc.z = 0;
    int idCell = (i * numeroColunasGridCell) + j;

    int kInicial = offSet[idCell];
    int kFinal = offSet[idCell] + count[idCell];

    //cout << "\t\t i  = " << i << "\t j = " << j << "\tidCell = " << idCell
    //     << "\t kInicial = " << kInicial << "\t kFinal = " << kFinal << endl;


    for (int k = kInicial; k < kFinal; k++) {
        int indiceParticulas = values[k];
        if (indiceParticulas != indiceAtual) {
            //cout << "\t\t Vizinho (visc.) => " << indiceParticulas << endl;
            float distancia = sqrt((posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x) * (posicoesParticulasT[indiceAtual].x - posicoesParticulasT[indiceParticulas].x) + (posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y) * (posicoesParticulasT[indiceAtual].y - posicoesParticulasT[indiceParticulas].y));
            if (distancia > 0 && distancia < h) {
                //cout << "\t\t (Distancia menor que h)" << endl;
                forcaVisc.x += massa * ((particulasT[indiceParticulas].velocidade.x + particulasT[indiceAtual].velocidade.x) / 2 * particulasT[indiceParticulas].massaEspecifica) * (45 / 3.14 * pow(h, 6) * h - distancia);
                forcaVisc.y += massa * ((particulasT[indiceParticulas].velocidade.y + particulasT[indiceAtual].velocidade.y) / 2 * particulasT[indiceParticulas].massaEspecifica) * (45 / 3.14 * pow(h, 6) * h - distancia);
                //forcaVisc.z += massa * ((particulasT[indiceParticulas].velocidade.z + particulasT[indiceAtual].velocidade.z) / 2 * particulasT[indiceParticulas].massaEspecifica) * (45/3.14*pow(h,6)* h - distanciaZ);
            }
        }
    }

    //cout << "\t\t Valor da viscosidade.x foi " << forcaVisc.x << "\t viscosidade.y foi " << forcaVisc.y << "\t viscosidade.z foi " << forcaVisc.z << endl;

    return forcaVisc;
}

void corrigePosicao(int posicaoParticulaAtual){
    
    if(posicoesParticulasT2[posicaoParticulaAtual].x < (XminGrid + recuoXGrid)){
        float novaPosicao = ABS((XminGrid + recuoXGrid) - posicoesParticulasT2[posicaoParticulaAtual].x);
        if(novaPosicao > ((XmaxGrid - recuoXGrid) - (XminGrid + recuoXGrid))){
            novaPosicao = (XmaxGrid - recuoXGrid) - (XminGrid + recuoXGrid);
        }
        //posicoesParticulasT2[posicaoParticulaAtual].x = (XminGrid + recuoXGrid) + novaPosicao;
        posicoesParticulasT2[posicaoParticulaAtual].x = (XmaxGrid - recuoXGrid) - novaPosicao;
        //posicoesParticulasT2[posicaoParticulaAtual].x = (XmaxGrid - recuoXGrid - recuoXGrid);
        
    }
    if(posicoesParticulasT2[posicaoParticulaAtual].x > (XmaxGrid - recuoXGrid)){
        float novaPosicao = ABS(posicoesParticulasT2[posicaoParticulaAtual].x - (XmaxGrid - recuoXGrid));
        if(novaPosicao > ((XmaxGrid - recuoXGrid) - (XminGrid + recuoXGrid))){
            novaPosicao = (XmaxGrid - recuoXGrid) - (XminGrid + recuoXGrid);
        }
        //posicoesParticulasT2[posicaoParticulaAtual].x = (XmaxGrid - recuoXGrid) - novaPosicao;
        posicoesParticulasT2[posicaoParticulaAtual].x = (XminGrid + recuoXGrid) + novaPosicao;
        //posicoesParticulasT2[posicaoParticulaAtual].x = (XminGrid + recuoXGrid + recuoXGrid);
    }
    if(posicoesParticulasT2[posicaoParticulaAtual].y < (YminGrid + recuoYGrid)){
        float novaPosicao = ABS((YminGrid + recuoYGrid) - posicoesParticulasT2[posicaoParticulaAtual].y);
        if(novaPosicao > ((YmaxGrid - recuoYGrid) - (YminGrid + recuoYGrid))){
            novaPosicao = (YmaxGrid - recuoYGrid) - (YminGrid + recuoYGrid);
        }
        //posicoesParticulasT2[posicaoParticulaAtual].y = (YminGrid + recuoYGrid) + novaPosicao;
        posicoesParticulasT2[posicaoParticulaAtual].y = (YmaxGrid - recuoYGrid) - novaPosicao;
        //posicoesParticulasT2[posicaoParticulaAtual].y = (YmaxGrid - recuoYGrid - recuoYGrid);
    }
    if(posicoesParticulasT2[posicaoParticulaAtual].y > (YmaxGrid - recuoYGrid)){
        float novaPosicao = ABS(posicoesParticulasT2[posicaoParticulaAtual].y - (YmaxGrid - recuoYGrid));
        if(novaPosicao > ((YmaxGrid - recuoYGrid) - (YminGrid + recuoYGrid))){
            novaPosicao = (YmaxGrid - recuoYGrid) - (YminGrid + recuoYGrid);
        }
        //posicoesParticulasT2[posicaoParticulaAtual].y = (YmaxGrid - recuoYGrid) - novaPosicao;
        posicoesParticulasT2[posicaoParticulaAtual].y = (YminGrid + recuoYGrid) + novaPosicao;
        //posicoesParticulasT2[posicaoParticulaAtual].y = (YminGrid + recuoYGrid + recuoYGrid);
    }

}

void calculaKernels(void) {

    if(VBO_Ativo){
        posicoesParticulasT = posicoesParticulas;
        posicoesParticulasT2 = posicoesParticulas2;
        particulasT = particulas;
        particulasT2 = particulas2;
    } else {
        posicoesParticulasT = posicoesParticulas2;
        posicoesParticulasT2 = posicoesParticulas;
        particulasT = particulas2;
        particulasT2 = particulas;
    }

    
    if (VBO_Ativo) {
        // Map buffer object for writing from CUDA
        CUDA_SAFE_CALL(cudaGLMapBufferObject((void**) &VBO, positionsVBO));

        launch_ZeraCountEOffSet(blocksPerGrid, threadsPerBlock, numeroCelulasGrid, countDevice, offSetDevice);

        launch_CalculaIDGridCell(blocksPerGrid, threadsPerBlock, VBO, valuesDevice,
                keysDevice, numeroParticulas, numeroColunasGridCell, deltaXGrid, deltaYGrid, countDevice,
                XminGrid, YminGrid);

        scan();

        // Copy result from device memory to host memory
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        //cudaMemcpy(values, valuesDevice, numeroParticulas*sizeof(unsigned int), cudaMemcpyDeviceToHost);
        CUDA_SAFE_CALL(cudaMemcpy(count, countDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(offSet, offSetDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));


#ifdef DEBUGCALC
        printCells();

        printCount();

        cout << "\n Fazendo o radixSort" << endl;
#endif

        radixSort();
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(values, valuesDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUGCALC
        for (int i = 0; i < numeroCelulasGrid; i++) {
            cout << "offSet: " << offSet[i] << "\t count[" << i << "] : " << count[i] << endl;
        }

        /*for (int i = 0; i < numeroParticulas; i++) {
            if (i < numeroCelulasGrid) {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << "\t offSet: " << offSet[i]
                        << "\t count[" << i << "] : " << count[i] << endl;
            } else {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << endl;
            }
        }*/

        printCells();
//        getchar();
#endif

        CUDA_SAFE_CALL(cudaMemcpy(particulas, particulasDevice, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(particulas2, particulasDevice2, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));

        // Unmap buffer object
        CUDA_SAFE_CALL(cudaGLUnmapBufferObject(positionsVBO));

    } else {

        // Map buffer object for writing from CUDA
        CUDA_SAFE_CALL(cudaGLMapBufferObject((void**) &VBO2, positionsVBO2));

        launch_ZeraCountEOffSet(blocksPerGrid, threadsPerBlock, numeroCelulasGrid, countDevice, offSetDevice);

        launch_CalculaIDGridCell(blocksPerGrid, threadsPerBlock, VBO2, valuesDevice,
                keysDevice, numeroParticulas, numeroColunasGridCell, deltaXGrid, deltaYGrid, countDevice,
                XminGrid, YminGrid);

        scan();

        // Copy result from device memory to host memory
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        //cudaMemcpy(values, valuesDevice, numeroParticulas*sizeof(unsigned int), cudaMemcpyDeviceToHost);
        CUDA_SAFE_CALL(cudaMemcpy(count, countDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(offSet, offSetDevice, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUGCALC
        printCells();

        printCount();

        cout << "\n Fazendo o radixSort" << endl;
#endif

        radixSort();
        CUDA_SAFE_CALL(cudaMemcpy(keys, keysDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(values, valuesDevice, numeroParticulas * sizeof (unsigned int), cudaMemcpyDeviceToHost));

#ifdef DEBUGCALC
        for (int i = 0; i < numeroCelulasGrid; i++) {
            cout << "offSet: " << offSet[i] << "\t count[" << i << "] : " << count[i] << endl;
        }

        /*for (int i = 0; i < numeroParticulas; i++) {
            if (i < numeroCelulasGrid) {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << "\t offSet: " << offSet[i]
                        << "\t count[" << i << "] : " << count[i] << endl;
            } else {
                cout << "id: " << values[i] << "\t idCell: " << keys[i] << endl;
            }
        }*/

        printCells();
//        getchar();
#endif

        CUDA_SAFE_CALL(cudaMemcpy(particulas, particulasDevice, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));
        CUDA_SAFE_CALL(cudaMemcpy(particulas2, particulasDevice2, numeroParticulas * sizeof (Particula), cudaMemcpyDeviceToHost));

        // Unmap buffer object
        CUDA_SAFE_CALL(cudaGLUnmapBufferObject(positionsVBO2));

    }

    deleteVBO();
    //glBindBufferARB(GL_ARRAY_BUFFER_ARB, positionsVBO);
    //float3* ptr = (float3*) glMapBufferARB(GL_ARRAY_BUFFER_ARB, GL_READ_WRITE_ARB);

    for (int indiceValue = 0; indiceValue < numeroParticulas; indiceValue++) {
        int jAtual = 0;
        int iAtual = 0;
        int jMin = 0;
        int jMax = numeroColunasGridCell;
        int iMin = 0;
        int iMax = numeroLinhasGridCell;
        //float h = deltaXGrid * 1.1; //sqrt(pow(deltaXGrid, 2) + pow(deltaYGrid, 2));
        float h = passoX/2 * 1.02;
        if (passoX > passoY) {
            h = passoY/2 * 1.02;
        }
        float3 forcaPressao = {0.0f, 0.0f, 0.0f};
        float3 forcaAceleracao = {0.0f, 0.0f, 0.0f};
        float3 forcaViscosidade = {0.0f, 0.0f, 0.0f};
        float3 gravidade = {0.0f, -9.8f, 0.0f};
        /*forcaPressaoVetor[indiceValue].x = forcaPressaoVetor[indiceValue].y = forcaPressaoVetor[indiceValue].z = 0;
        forcaAceleracaoVetor[indiceValue].x = forcaAceleracaoVetor[indiceValue].y = forcaAceleracaoVetor[indiceValue].z = 0;
        forcaViscosidadeVetor[indiceValue].x = forcaViscosidadeVetor[indiceValue].y = forcaViscosidadeVetor[indiceValue].z = 0;*/

        cout << "\nPara o indice indice value igual a = " << indiceValue << " e values[indiceValue] = " << values[indiceValue] << endl;

        if (indiceValue < numeroParticulas) {
            int posicaoParticulaAtual = values[indiceValue];
            jAtual = (posicoesParticulasT[posicaoParticulaAtual].x / deltaXGrid);
            iAtual = (posicoesParticulasT[posicaoParticulaAtual].y / deltaYGrid);
            //int idCellAtual = (iAtual * numeroColunasGridCell) + jAtual;
            cout << "\t jAtual = " << jAtual << "    iAtual = " << iAtual;

            //VERIFICAO DAS EXTREMIDADES.
            if (jAtual == 0) {
                jMin = 0;
            } else {
                jMin = jAtual - 1;
            }

            if (jAtual == (numeroColunasGridCell - 1)) {
                jMax = jAtual;
            } else {
                jMax = jAtual + 1;
            }

            if (iAtual == 0) {
                iMin = 0;
            } else {
                iMin = iAtual - 1;
            }

            if (iAtual == (numeroLinhasGridCell - 1)) {
                iMax = iAtual;
            } else {
                iMax = iAtual + 1;
            }

            cout << "\t (jMin = " << jMin << "    jMax = " << jMax << "    iMin = " << iMin << "    iMax = " << iMax << ")" << endl;
            particulasT2[posicaoParticulaAtual].massaEspecifica = 0;

            //CALCULA AS MASSAS ESPECIFICAS.
            for (int i = iMin; i <= iMax; i++) {
                for (int j = jMin; j <= jMax; j++) {
                    particulasT2[posicaoParticulaAtual].massaEspecifica += calculaMassaEspecifica(i, j, h, posicaoParticulaAtual);
                    //cout << "\t\t Massa esp. temp da particula " << posicaoParticulaAtual << " eh " << particulasT1[posicaoParticulaAtual].massaEspecifica << endl;
                }
            }

            cout << "\t pos.x: " << posicoesParticulasT[posicaoParticulaAtual].x
                    << "\t pos.y: " << posicoesParticulasT[posicaoParticulaAtual].y
                    << "\t pos.z: " << posicoesParticulasT[posicaoParticulaAtual].z << endl;
            cout << "\t Massa especifica final da particula " << posicaoParticulaAtual << " eh " << particulasT2[posicaoParticulaAtual].massaEspecifica << endl ;
            
            //CALCALA A PRESSAO NA PARTICULA.
            particulasT2[posicaoParticulaAtual].pressao = constGas * (particulasT2[posicaoParticulaAtual].massaEspecifica - massaEspecificaInicial);

            //CALCULA A FORÇA DE PRESSAO.
            for (int i = iMin; i <= iMax; i++) {
                for (int j = jMin; j <= jMax; j++) {
                    float3 forcaPressaoTemp = calculaFocaPressao(i, j, h, posicaoParticulaAtual);
                    forcaPressao.x += forcaPressaoTemp.x;
                    forcaPressao.y += forcaPressaoTemp.y;
                    forcaPressao.z += forcaPressaoTemp.z;
                }
            }
            forcaPressao.x *= -1;
            forcaPressao.y *= -1;
            forcaPressao.z *= -1;
            forcaPressaoVetor[posicaoParticulaAtual].x = forcaPressao.x;
            forcaPressaoVetor[posicaoParticulaAtual].y = forcaPressao.y;
            forcaPressaoVetor[posicaoParticulaAtual].z = forcaPressao.z;
            cout << "\t ForcaPressao.x: " << forcaPressao.x << "\t ForcaPressao.y: " << forcaPressao.y << "\t ForcaPressao.z: " << forcaPressao.z << endl;

            //CALCULA A FORÇA VISCOSIDADE.
            for (int i = iMin; i <= iMax; i++) {
                for (int j = jMin; j <= jMax; j++) {
                    float3 forcaViscTemp = calculaFocaVisc(i, j, h, posicaoParticulaAtual);
                    forcaViscosidade.x += forcaViscTemp.x;
                    forcaViscosidade.y += forcaViscTemp.y;
                    forcaViscosidade.z += forcaViscTemp.z;
                }
            }
            forcaViscosidade.x *= coeficienteViscosidade;
            forcaViscosidade.y *= coeficienteViscosidade;
            forcaViscosidade.z *= coeficienteViscosidade;
            forcaViscosidadeVetor[posicaoParticulaAtual].x = forcaViscosidade.x;
            forcaViscosidadeVetor[posicaoParticulaAtual].y = forcaViscosidade.y;
            forcaViscosidadeVetor[posicaoParticulaAtual].z = forcaViscosidade.z;
            cout << "\t ForcaViscosidade.x: " << forcaViscosidade.x << "\t ForcaViscosidade.y: " << forcaViscosidade.y << "\t ForcaViscosidade.z: " << forcaViscosidade.z << endl;
            
            forcaAceleracao.x = forcaPressao.x + forcaViscosidade.x; //+ particulasT1[posicaoParticulaAtual].massaEspecifica * gravidade.x;
            forcaAceleracao.y = forcaPressao.y + forcaViscosidade.y; //+ particulasT1[posicaoParticulaAtual].massaEspecifica * gravidade.y;
            forcaAceleracao.z = forcaPressao.z + forcaViscosidade.z; //+ particulasT1[posicaoParticulaAtual].massaEspecifica * gravidade.z;
            forcaAceleracaoVetor[posicaoParticulaAtual].x = forcaAceleracao.x;
            forcaAceleracaoVetor[posicaoParticulaAtual].y = forcaAceleracao.y;
            forcaAceleracaoVetor[posicaoParticulaAtual].z = forcaAceleracao.z;
            cout << "\t ForcaAceleracao.x: " << forcaAceleracao.x << "\t ForcaAceleracao.y: " << forcaAceleracao.y << "\t ForcaAceleracao.z: " << forcaAceleracao.z << endl;

            particulasT2[posicaoParticulaAtual].velocidade.x = particulasT[posicaoParticulaAtual].velocidade.x + passoIntegracao / 2 * forcaAceleracao.x;
            particulasT2[posicaoParticulaAtual].velocidade.y = particulasT[posicaoParticulaAtual].velocidade.y + passoIntegracao / 2 * forcaAceleracao.y;
            particulasT2[posicaoParticulaAtual].velocidade.z = 0;//particulasT[posicaoParticulaAtual].velocidade.z + passoIntegracao / 2 * forcaAceleracao.z;
            cout << "\t Velocidade.x: " << particulasT2[posicaoParticulaAtual].velocidade.x
                    << "\t Velocidade.y: " << particulasT2[posicaoParticulaAtual].velocidade.y
                    << "\t Velocidade.z: " << particulasT2[posicaoParticulaAtual].velocidade.z << endl;

            posicoesParticulasT2[posicaoParticulaAtual].x = posicoesParticulasT[posicaoParticulaAtual].x + passoIntegracao / 2 * particulasT2[posicaoParticulaAtual].velocidade.x;
            cout << passoIntegracao / 2 * particulasT2[posicaoParticulaAtual].velocidade.x;
            posicoesParticulasT2[posicaoParticulaAtual].y = posicoesParticulasT[posicaoParticulaAtual].y + passoIntegracao / 2 * particulasT2[posicaoParticulaAtual].velocidade.y;
            cout << "\t" << passoIntegracao / 2 * particulasT2[posicaoParticulaAtual].velocidade.x << endl;
            posicoesParticulasT2[posicaoParticulaAtual].z = 0;//posicoesParticulasT[posicaoParticulaAtual].z + passoIntegracao / 2 * particulasT1[posicaoParticulaAtual].velocidade.z;

            corrigePosicao(posicaoParticulaAtual);

            cout << "\t pos.x: " << posicoesParticulasT2[posicaoParticulaAtual].x
                    << "\t pos.y: " << posicoesParticulasT2[posicaoParticulaAtual].y
                    << "\t pos.z: " << posicoesParticulasT2[posicaoParticulaAtual].z << endl;

            /*int jTemp = (posicoesParticulasT2[posicaoParticulaAtual].x) / deltaXGrid;
            int iTemp = (posicoesParticulasT2[posicaoParticulaAtual].y) / deltaYGrid;
            int idCell = (iTemp * numeroColunasGridCell) + jTemp;
            cout << indiceValue << "-> Next idCell: " << idCell << endl;*/

        }
    }

    if (VBO_Ativo) {
        VBO_Ativo = false;
    } else {
        VBO_Ativo = true;
    }

    buildVBO();

    glutPostRedisplay();
}

void teclado(unsigned char key, int x, int y) {

    switch (key) {

        case 27:
            exit(0);
            break;
        case 'g':
        case 'G':
            if (mostraGrid) {
                mostraGrid = false;
#ifdef DEBUG
                cout << "MostraGrid = false;" << endl;
#endif
            } else {
                mostraGrid = true;
#ifdef DEBUG
                cout << "MostraGrid = true;" << endl;
#endif

            }
            break;
        case 'e':
        case 'E':
            if (mostraEixo) {
                mostraEixo = false;
#ifdef DEBUG
                cout << "MostraEixo = false;" << endl;
#endif

            } else {
                mostraEixo = true;
#ifdef DEBUG
                cout << "MostraEixo = true;" << endl;
#endif

            }
            break;
        case 'm':
        case 'M':
            if (mostraVetores) {
                mostraVetores = false;
#ifdef DEBUG
                cout << "MostraVetores = false;" << endl;
#endif

            } else {
                mostraVetores = true;
#ifdef DEBUG
                cout << "MostraVetores = true;" << endl;
#endif

            }
            break;
        case 'n':
        case 'N':
            if (mostraPressao) {
                mostraPressao = false;
#ifdef DEBUG
                cout << "MostraPressao = false;" << endl;
#endif

            } else {
                mostraPressao = true;
#ifdef DEBUG
                cout << "MostraPressao = true;" << endl;
#endif

            }
            break;
        case 'b':
        case 'B':
            if (mostraViscosidade) {
                mostraViscosidade = false;
#ifdef DEBUG
                cout << "MostraViscosidade = false;" << endl;
#endif

            } else {
                mostraViscosidade = true;
#ifdef DEBUG
                cout << "MostraViscosidade = true;" << endl;
#endif

            }
            break;
            case 'v':
            case 'V':
                mostraAceleracao = !mostraAceleracao;
#ifdef DEBUG
                cout << "MostraViscosidade = " << mostraAceleracao << endl;
#endif

                break;
        case 's':
        case 'S':
            if (gpu) {
                CalculaKernelsOnGPU();
            } else {
                calculaKernels();
                /*if (VBO_Ativo) {
                    cout << "Imprimindo particulas t (particulas)." << endl;
                    printDadosParticulas();
                    cout << "\n Imprimindo particulas t + 1 (particulas2)." << endl;
                    printDadosParticulas2();
                } else {
                    cout << "Imprimindo particulas t + 1 (particulas)." << endl;
                    printDadosParticulas();
                    cout << "\n Imprimindo particulas t (particulas2)." << endl;
                    printDadosParticulas2();
                }*/
            }
            break;
        case 'r':
        case 'R':
            cout << "Processo de simulação estartado!" << endl;
            glutIdleFunc(calculaKernels);
            break;
        case 'p':
        case 'P':
            glutIdleFunc(NULL);
            break;
    }
            
    glutPostRedisplay();

}

void iniciaDados() {

    mostraGrid = true;
    mostraEixo = false;
    mostraVetores = true;
    mostraPressao = true;
    mostraViscosidade = true;
    mostraAceleracao = true;

    //deltaXGrid = (float *) malloc(sizeof (float));
    //deltaYGrid = (float *) malloc(sizeof (float));

    deltaXGrid = (XmaxGrid - XminGrid) / numeroColunasGridCell;
    deltaYGrid = (YmaxGrid - YminGrid) / numeroLinhasGridCell;

    numeroCelulasGrid = numeroColunasGridCell * numeroLinhasGridCell;
    count = (unsigned int *) malloc(numeroCelulasGrid * sizeof (unsigned int));
    offSet = (unsigned int *) malloc(numeroCelulasGrid * sizeof (unsigned int));
    //Aloca o vetor de posições na memoria do host
    posicoesParticulas = (float3 *) malloc(numeroParticulas * sizeof (float3));
    posicoesParticulas2 = (float3 *) malloc(numeroParticulas * sizeof (float3));
    forcaPressaoVetor = (float3 *) malloc(numeroParticulas * sizeof (float3));
    forcaAceleracaoVetor = (float3 *) malloc(numeroParticulas * sizeof (float3));
    forcaViscosidadeVetor = (float3 *) malloc(numeroParticulas * sizeof (float3));
    particulas = (Particula *) malloc(numeroParticulas * sizeof (Particula));
    particulas2 = (Particula *) malloc(numeroParticulas * sizeof (Particula));
    keys = (unsigned int *) malloc(numeroParticulas * sizeof (unsigned int));
    values = (unsigned int *) malloc(numeroParticulas * sizeof (unsigned int));

    linhas = numeroParticulas / colunas;

    passoX = ((XmaxGrid - recuoXGrid) - (XminGrid + recuoXGrid)) / colunas;
    passoY = ((YmaxGrid - recuoYGrid) - (YminGrid + recuoYGrid)) / linhas;

    posicaoX = XminGrid + recuoXGrid + 0.23;
    posicaoY = YmaxGrid - recuoYGrid - 0.23;

    //Inicia os dados das particulas e gera os dados.
    for (int i = 0; i < numeroParticulas; i++) {
        particulas[i].massaEspecifica = massaEspecificaInicial;
        particulas[i].pressao = 1;
        particulas[i].velocidade.x = 0;
        particulas[i].velocidade.y = 0;
        particulas[i].velocidade.z = 0;
        particulas2[i].massaEspecifica = massaEspecificaInicial;
        particulas2[i].pressao = 1;
        particulas2[i].velocidade.x = 0;
        particulas2[i].velocidade.y = 0;
        particulas2[i].velocidade.z = 0;
        //Gera os values.
        values[i] = i;
        forcaPressaoVetor[i].x = forcaPressaoVetor[i].y = forcaPressaoVetor[i].z = 0;
        forcaAceleracaoVetor[i].x = forcaAceleracaoVetor[i].y = forcaAceleracaoVetor[i].z = 0;
        forcaViscosidadeVetor[i].x = forcaViscosidadeVetor[i].y = forcaViscosidadeVetor[i].z = 0;
    }

    threadsPerBlock = 512;
    blocksPerGrid = (numeroParticulas + threadsPerBlock - 1) / threadsPerBlock;
}

void printDadosIniciais() {
#ifdef DEBUG
    cout << "\n-------------------- DADOS DA SIMULACAO --------------------" << endl;
    cout << "Numero de particulas: " << numeroParticulas << endl;
    cout << "Numero de linhas do Grid: " << numeroLinhasGridCell << endl;
    cout << "Numero de colunas do Grid: " << numeroColunasGridCell << endl;
    cout << "Numero de colunas para distribuicao das particulas: " << colunas << endl;
    cout << "DeltaX: " << deltaXGrid << endl;
    cout << "DeltaY: " << deltaYGrid << endl;
    cout << "------------------------------------------------------------" << endl;
//    getchar();
#endif
}

void iniciaParticulas() {

    //PREENCHENDO AS POSIÇÕES X.
    int contElementosLinha = 0;
    for (int j = 0; j < numeroParticulas; j++) {
        posicoesParticulas[j].x = posicaoX;
        posicoesParticulas2[j].x = posicaoX;

        contElementosLinha++;
        posicaoX = posicaoX + passoX/2;
        if (contElementosLinha >= colunas) {
            posicaoX = XminGrid + recuoXGrid + 0.23;
            contElementosLinha = 0;
        }
    }

    //PREENCHENDO AS POSIÇÕES Y.
    int contElementosColuna = 0;
    for (int j = 0; j < numeroParticulas; j++) {
        posicoesParticulas[j].y = posicaoY;
        posicoesParticulas2[j].y = posicaoY;
        contElementosColuna++;
        if (contElementosColuna >= colunas) {
            posicaoY = posicaoY - passoY/2;
            contElementosColuna = 0;
        }
    }

    //PREENCHENDO AS POSIÇÕES Z.
    for (int j = 0; j < numeroParticulas; j++) {
            posicoesParticulas[j].z = 0.0;
            posicoesParticulas2[j].z = 0.0;
    }
}

void iniciaCuda() {

    //Alocando as constantes da aplicação.
    /*CUDA_SAFE_CALL(cudaMalloc((void **) &massa_d, sizeof (float)));
    CUDA_SAFE_CALL(cudaMemcpyToSymbol("massa_d", &massa, sizeof (float), 0, cudaMemcpyHostToDevice));
    CUDA_SAFE_CALL(cudaGetSymbolAddress((void **)&massa_d, "massa_d"));*/
    //Aloca o vetor de particulas na memoria do device
    CUDA_SAFE_CALL(cudaMalloc((void **) & particulasDevice, numeroParticulas * sizeof (Particula)));
    CUDA_SAFE_CALL(cudaMalloc((void **) & particulasDevice2, numeroParticulas * sizeof (Particula)));
    CUDA_SAFE_CALL(cudaMalloc((void **) & keysDevice, numeroParticulas * sizeof (unsigned int)));
    CUDA_SAFE_CALL(cudaMalloc((void **) & valuesDevice, numeroParticulas * sizeof (unsigned int)));
    CUDA_SAFE_CALL(cudaMalloc((void **) & countDevice, numeroCelulasGrid * sizeof (unsigned int)));
    CUDA_SAFE_CALL(cudaMalloc((void **) & offSetDevice, numeroCelulasGrid * sizeof (unsigned int)));
    //cudaMalloc((void **) & deltaXGridDevice, sizeof (float));
    //cudaMalloc((void **) & deltaYGridDevice, sizeof (float));

    CUDA_SAFE_CALL(cudaMemcpy(valuesDevice, values, numeroParticulas * sizeof (unsigned int), cudaMemcpyHostToDevice));
    CUDA_SAFE_CALL(cudaMemcpy(particulasDevice, particulas, numeroParticulas * sizeof (Particula), cudaMemcpyHostToDevice));
    CUDA_SAFE_CALL(cudaMemcpy(particulasDevice2, particulas2, numeroParticulas * sizeof (Particula), cudaMemcpyHostToDevice));
    //cudaMemcpy(keysDevice, keys, numeroParticulas*sizeof(unsigned int), cudaMemcpyHostToDevice);
    CUDA_SAFE_CALL(cudaMemcpy(countDevice, count, numeroCelulasGrid * sizeof (unsigned int), cudaMemcpyHostToDevice));
    //cudaMemcpy(deltaXGridDevice, deltaXGrid, sizeof (float), cudaMemcpyHostToDevice);
    //cudaMemcpy(deltaYGridDevice, deltaYGrid, sizeof (float), cudaMemcpyHostToDevice);


}

/*
 *
 */
int main(int argc, char** argv) {

    //Verifica se todos os parametros realmente foram passados na chamada do programa.
    if (argc < 5) {
        cout << "Falta parametros!" << endl;
        exit(EXIT_SUCCESS);
    }

    numeroParticulas = atof(argv[1]);
    numeroColunasGridCell = atof(argv[2]);
    numeroLinhasGridCell = atof(argv[3]);
    colunas = atof(argv[4]);
    massa = 50/numeroParticulas;
    cout << "Massa = " << massa << "g." << endl;
    massaEspecificaInicial = (0.5/0.25);
    cout << "Massa especifica inicial = " << massaEspecificaInicial << "g." << endl;

    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
    glutInitWindowSize(800, 600);
    glutInitWindowPosition(10, 10);
    glutCreateWindow("Simulador de fluido 2D");
    GLenum err = glewInit();
    if (GLEW_OK != err) {
        //Problem: glewInit failed, something is seriously wrong.
        fprintf(stderr, "Error: %s\n", glewGetErrorString(err));
    }

    //printDeviceProps();
    cudaGLSetGLDevice(0);

    iniciaDados();
#ifdef DEBUG
    printDadosIniciais();
    printDadosParticulas();
#endif

    iniciaParticulas();

#ifdef DEBUG
    printPosicoesParticulas();
#endif

#ifdef DEBUG_GL
    cout << "ANTES DO VBO!" << endl;
#endif
    buildVBO();
    VBO_Ativo = true;
#ifdef DEBUG_GL
    cout << "DEPOIS DO VBO!" << endl;
#endif

    posicoesParticulasT = posicoesParticulas;
    posicoesParticulasT2 = posicoesParticulas2;
    particulasT = particulas;
    particulasT2 = particulas2;

    init();
    iniciaCuda();

    //OPENGL
    glutKeyboardFunc(teclado);
    glutDisplayFunc(display);
    glutReshapeFunc(AlteraTamanhoJanela);
    //glutIdleFunc(Desenha);

    glutMainLoop();

    return (EXIT_SUCCESS);
}
